Cumulative human impacts on global marine fauna highlight risk to biological and functional diversity

Anthropogenic stressors to marine ecosystems from climate change and human activities increase extinction risk of species, disrupt ecosystem integrity, and threaten important ecosystem services. Addressing these stressors requires understanding where and to what extent they are impacting marine biological and functional diversity. We model cumulative risk of human impact upon 21,159 marine animal species by combining information on species-level vulnerability and spatial exposure to a range of anthropogenic stressors. We apply this species-level assessment of human impacts to examine patterns of species-stressor interactions within taxonomic groups. We then spatially map impacts across the global ocean, identifying locations where climate-driven impacts overlap with fishing, shipping, and land-based stressors to help inform conservation needs and opportunities. Comparing species-level modeled impacts to those based on marine habitats that represent important marine ecosystems, we find that even relatively untouched habitats may still be home to species at elevated risk, and that many species-rich coastal regions may be at greater risk than indicated from habitat-based methods alone. Finally, we incorporate a trait-based metric of functional diversity to identify where impacts to functionally unique species might pose greater risk to community structure and ecosystem integrity. These complementary lenses of species, function, and habitat provide a richer understanding of threats to marine biodiversity to help inform efforts to meet conservation targets and ensure sustainability of nature’s contributions to people.


Supporting Figures
. Graphical conceptual summary of cumulative impact calculations.(A) Stressors are mapped as intensity values on a 10 x 10 km grid, with values ranging from 0 to 1. (B) Each species is vulnerable, to different degrees, to each of the stressors, with vulnerability ranging from 0 to 1. (C) The distribution of each species is mapped as present/absent in each 10 x 10 km grid cell.Note the colors represent different functional entities.(D) For each species x stressor combination, an impact map is generated as the product of species presence (C), species vulnerability to the stressor (B), and stressor intensity (A).(E) The cumulative impact to each species is the sum of impacts across all stressors.(F) For the species approach, a mean cumulative human impact is calculated as the average cumulative impact across all species present in each cell.(G) For the functional entity approach, impact values are aggregated to the functional entity level by averaging impact values within each grid cell across all species associated with that functional entity.A functional vulnerability is calculated for each FE/grid cell combination, where a FE with few species is more vulnerable than a FE represented by many species.(H) For the functional entity approach, a mean cumulative human impact is calculated as the average cumulative impact across all functional entities present in each cell, weighted by functional vulnerability.(I) The functional entity impacts and species impacts are compared by taking the difference between the impacts predicted by each method for each grid cell.Comparison of cumulative, climate, and non-climate stressors by habitat and species methods across 10 km resolution cells within coastal ( !,  ! ) and oceanic ( " ,  " ) portions of 62 representative marine ecological provinces, transformed to percentile ranks relative to global distribution within each impact category.Filled point indicates median value; bars represent interquartile range (IQR, quartile Q1 to Q3); whiskers indicate observations 1.5x IQR below (above) Q1 (Q3) of bar.Outliers omitted from plots for clarity.A) Cumulative impacts by species and habitat cumulative impact methods.B) Climate impacts by species and habitat cumulative impact methods.C) Non-climate impacts by species and habitat cumulative impact methods.

Figure S6
. Sensitivity of functional vulnerability (FV) and cumulative human impact (CHI) to trait variation.Each trait was resampled 1000 times, holding other traits constant, and change in functional vulnerability (A-D) and cumulative human impact (E-H) were recalculated according to the Functional Entity approach.Functional vulnerability generally increased slightly (A-D, green tones), as resampling tends to create novel trait combinations not found in nature, leading to more low-membership (and thus high FV) functional entities.The effects of these changes in functional vulnerability on the overall cumulative impact were quite low.

Figure S7.
Trait imputation.A) Missingness of traits for inclusion in functional entities, prior to MICE imputation.Water column position (wcol), log(body length) (log_l), trophic level (troph), and adult mobility (adult_mob) were used to assign species to functional entities; dark boxes represent known values, yellow boxes represent missing values.Two other traits highlighted in the red box, log(fecundity) (log_f) and age to maturity (age_mat) were used to assist in imputation where available; color scale represents proportional missingness of those traits for each combination of the other four.B) After MICE imputation, an additional gapfill step using taxonomic relatives was used."None" represents proportion of species successfully imputed by MICE; "genus" represents proportion of species whose traits were gapfilled using the most common values of other species in its genus; and "family" represents the proportion filled using the most common values of other species in its family.

Analysis grid
All spatial analyses were calculated on a gridded global map using a Mollweide equal-area projection coordinate reference system (CRS), gridded to 10 km x 10 km resolution.An ocean base map was prepared by rasterizing the vector ocean polygon features of the Natural Earth (https://www.naturalearthdata.com/)Oceans 1:10m dataset to a 1 km x 1 km Mollweide projection, then aggregating by a factor of 10 to approximate percentage of ocean within each cell.The resulting 10 km Mollweide ocean raster was used as the target raster for projecting all other datasets, and was used to mask out non-ocean cells from reprojected data.
To examine variation in calculations across coastal regions and ecoregions, we prepared two supporting rasters at the same CRS as our analysis grid: coastal area based on cells with a minimum depth less than 200 m, based on GEBCO bathymetry data [90], and marine ecoregions of the world [46].

Species distributions
Species distribution data were taken from AquaMaps [30] (n = 18,480) and IUCN species distribution maps [31,32] (n = 2,679).For both datasets, synonymous scientific binomials and slight differences in nomenclature were resolved by comparing names against accepted names of marine species in the World Register of Marine Species (WoRMS, [29]) using the taxize package [43,44].For species appearing in both distribution map datasets, the AquaMaps distribution maps, based on transparent and repeatable algorithms using publicly available data, were preferred over IUCN range maps, which integrate data and expert knowledge but may include mapping decisions that are difficult to replicate.Each mapping method has advantages and disadvantages, though neither method can perfectly replicate the underlying truth [91].
AquaMaps predicts species ranges as "probability of occurrence" in 0.5° cells, computed based on long-term averages of species-specific relative environmental suitability across multiple environmental variables (temperature, salinity, depth, sea ice concentration, primary production, and in some cases oxygen level and distance to shore) based on species-specific habitat preference derived from publicly available occurrence data [30].To determine the range map for each species, we converted gradient probabilities into binary presence/absence using a threshold of ≥0.5 to represent "presence" of a given species, following a common practice in conservation literature [e.g., 50, [92][93][94].Species maps based on fewer than 10 "occurrence cells" (0.5° cells with at least one validated occurence record) were rejected to ensure quality [30].These results were then reprojected in raster form to the Mollweide CRS at 10 km resolution, resulting in n = 18,480 species presence maps.AquaMaps outputs have been validated against independent survey occurrence data and other algorithms often used in species distribution modeling [95,96].
The IUCN dataset presents species ranges as polygons representing the historical, present and possible distribution of a taxon's occurrences [31].For each species (or subpopulation where available), we excluded polygons with a "presence" value of 5 indicating "extinct" portions of a range, reprojected the remaining polygon features to the Mollweide CRS, then rasterized the results to the 10 km analysis grid, resulting in n = 2,679 species presence maps.
For coastal and neritic species from both datasets, we masked the resulting presence maps to cells with a minimum depth of 200 m or less using GEBCO bathymetry data [90], and masked using the ocean area raster to exclude non-ocean cells.The depth mask for shallow-water species is an effective method of reducing errors of commission in IUCN range maps to better align them with AquaMaps distribution maps [97].While AquaMaps uses bottom depth as a parameter to limit ranges for depth-limited species, this masking method allows for a finer resolution map than the 0.5 degree native resolution of AquaMaps.

Vulnerability estimates
Vulnerability weights, i.e., relative effect of stressor  on the fitness/health of species , were determined based on methods of Butt et al. [19].That study estimated vulnerability of species  to stressor  based on presence of certain traits that are likely to increase the species' physiological sensitivity to the stressor  !" , other traits that affect the species' ability to adapt to that specific stressor, i.e., stressor-specific adaptive capacity  !" , and life history and population-level traits that affect the population's ability to adapt to or recover from disturbances in general, i.e., general adaptive capacity  ! .Traits were binned into nominal or ordinal categories and certain trait values corresponded to different scores for different stressors (see S3 Table ).An additional exposure modifier was included to account for possibility of exposure  !" ∈ {0,1} of species  to stressor , e.g., a mesopelagic species (depth below 200 m) will not be exposed to ship strikes, so  !" = 0.These metrics were combined to produce a vulnerability score  !" ∈ [0,1]: For the present study, we updated some aspects of methodology to improve imputation of species vulnerability for species with partial trait sets.First, for species whose traits were provided by taxon experts at genus level or higher rank in the trait set assembled by Butt et al. [19], certain traits (i.e., body length, fecundity, generation time, temperature tolerances, and depth preferences) were filled using species-level data (where available) from FishBase/SealifeBase [36,37] and AquaMaps [30], while extent of occurrence was determined from the species' distribution map as described above.
Remaining traits were taken as provided by the taxon experts for the original study.Second, rather than calculate sensitivity, adaptive capacity, and vulnerability scores for species groups then impute missing values based on the distribution of those scores (as done in the original study), we first imputed missing trait values based on frequency within the species' taxonomic neighbors, then scored the sensitivity, adaptive capacity, and vulnerability from the combination of species-level traits (from FishBase, SeaLifeBase, and AquaMaps) and those imputed traits.Note that these methodological updates resulted in small differences in the vulnerability scores between the published study and our present study: mean/sd difference of .0002± .005, with 98% of differences between -0.00009 and 0.00024.
The expert-elicited trait values used to calculate vulnerability were given as ordinal or nominal categorical values, generally with a single value attributed to each trait for species.While a single value does not account for large-scale interspecific variability or regional variations in traits, most traits are broadly applicable across a species, genus, or family (e.g., respiration structures of gills vs. lungs vs. diffusion), and for traits representing continuous values (e.g., body size, fecundity, age to maturity), the bins typically increase as orders of magnitude, wide enough to accommodate substantial variation without affecting the vulnerability estimate.
The traits used to calculate sensitivity, adaptive capacity, and exposure components of vulnerability for each stressor are included in S3 Table .The updated trait-based vulnerability methods and results, including the full matrix translating trait values to vulnerability components and the full dataset of species traits, can be found at https://github.com/mapping-marine-sppvuln/spp_vuln_framework.

Uniform exposure stressors
For most of the included stressors, exposure does not depend on species identity (though vulnerability to the stressor certainly might) and therefore we modeled exposure as uniform across all species.For these stressors, a single map of relative stressor intensity was created from gridded data using the following general process: • We reprojected raw intensity to Mollweide CRS at 10 km resolution • For stressors where marginal impact is expected to be decreasing with intensity (e.g., the hundredth hour of trawling in an area likely overlaps habitat already destroyed by the first hour of trawling), we applied a log transformation to the raw data.• For stressors whose distribution contains a small number of extreme outliers, we identified a reference point based on the 99.9th percentile; otherwise we assigned a reference point based on the maximum observed value.• Finally, we rescaled the data using the reference point to result in a distribution of stressor intensity ranging between zero and one.
Information on the data source, transformation, and reference point used for each stressor layer can be found in S4 Table.

Non-uniform exposure stressors
For several stressors, exposure intensity depends on species-specific information.These stressors include bycatch, targeted fishing, and sea surface temperature rise.

Bycatch stressor
The degree to which species are exposed to bycatch is dependent on their position in the water column.Pelagic species are unlikely to be swept up in a bottom trawl, while demersal species are unlikely to be swept up in a midwater trawl or purse seine.We prepared three bycatch layers, summing industrial and nonindustrial discards based on gear type listed by Watson [33] and Watson et al. [85]: • Benthic bycatch (affects species identified as benthic), based on gear types: trawl, dredge, and trap • Pelagic bycatch (affects species identified as pelagic), based on gear types: line (tuna and nontuna), longline (tuna and non-tuna), midwater trawl, seine, purse seine (tuna and non-tuna), gillnet, other • Both (affects species identified as benthopelagic or reef-associated): the average of benthic and pelagic bycatch layers.
Catch estimates in Watson [33] and Watson et al. [85] provide data on discards (industrial/nonindustrial) by gear type, presented in 0.5° cells.The discard values were summed across benthic or pelagic gear types, then the totals were normalized by cell ocean area resulting in an intensity of discarded catch, i.e., tonnes of bycatch per square kilometer.These intensity rasters were reprojected to the 10 km Mollweide CRS analysis grid.The intensity rasters were then adjusted by dividing by ln(NPP) according to water column position, to indicate that bycatch in a high productivity area is less problematic than the same amount of bycatch in a low productivity area.Surface NPP data were taken from Bio-ORACLE [86,87], mean sea surface net primary productivity (NPP) of carbon, g/m 3 /day.
Benthic NPP sums productivity at bottom depth and export flux (e.g., "marine snow") from the surface to bottom depth.Bottom NPP data were taken from Bio-ORACLE [86,87], mean NPP of carbon at mean bottom depth, g/m 3 /day.Export flux from surface to bottom depth were calculated based on an exponential decay model for export flux at depth : () =  ' × (1 − ) ( # Applying non-linear least squares using data from Table 1 (control) in Gt C a -1 (globally integrated) in Yool et al [98], we identified best fit parameters  = 0.341,  = .288.
Finally, the resulting surface and benthic catch/ln(NPP) layers were rescaled from 0 to 1, using a reference point based on the 99.9% quantile of observed cell values.

Targeted fishing stressor
In addition to discards, Watson [33] and Watson et al. [85] report targeted catch for industrial and non-industrial fisheries at 0.5° cells, across multiple gear types and taxonomic groups.While any species might be vulnerable to targeted fishing, not all species are targeted, and so the targeted fishing stressor layer is distinct for every targeted species (and nonexistent for non-targeted species), thus a targeted fishing stressor layer was calculated separately for each species with non-zero catch in the Watson dataset.
Taxon names were compared to accepted names per WoRMS [29] using the taxize package [43,44] to resolve synonyms and differences in spelling.Total catch for each taxon was summed across pelagic gears and benthic gears separately, then divided by cell ocean area resulting in intensity, i.e., tonnes of catch per km 2 .Catch reported at the species level was attributed directly to that species.Catch in a given cell but reported at higher ranks (e.g., genus, family) was divided equally among all local species (per species distributions) in that genus or family.In many cases, a given species would be attributed catch at multiple taxonomic levels in the same cell, which were summed to create a cell total catch intensity for each species, separately for pelagic and benthic gear types.The pelagic and benthic catch intensities were reprojected to the 10 km Mollweide CRS analysis grid, and then normalized by ln(NPP), either surface or benthic as appropriate, to account for the fact that a unit of catch in a highly productive region imposes less stress on an ecosystem than the same unit of catch in a low-productivity region.The pelagic and benthic NPP-normalized catch were then summed for each cell for the species.
Reference points to rescale the targeted fishing stressor layers are species specific.A global maximum reference point was set by first calculating the 90th percentile of NPP-normalized catch for each species across its entire distribution, then selecting the score of the species with the highest 90th percentile value: Engraulis ringens, Peruvian anchoveta, at  )*+ = 2,170 tonnes of NPPnormalized catch.This global reference point was used to rescale (from 0 to 1) any species whose 99.9th percentile of NPP-normalized catch exceeded this value (15 species total).For species whose 99.9th percentile catch across its range falls below this reference point, we used the 99.9th percentile of that species' catch across its range as its own reference point.The reference catch for species  is therefore: The NPP-normalized catch was then rescaled using the appropriate reference point, with values capped at 1.0, resulting in a gridded map of stressor values from 0 to 1 for every species with nonzero targeted catch.

SST rise stressor
We included two stressors related to ocean temperature: sea surface temperature extremes, representing impacts from short-lived (weeks to months) high temperature events, i.e., marine heat waves; and rise in annual mean sea surface temperature representing long-term (years to decades) changes in sea surface temperature relative to a historic baseline.The SST extremes stressor is described above in the uniform-exposure stressors.Exposure to long-term SST rise estimates the impact to a species when mean annual temperatures risk exceeding the physiological tolerance of the species due to the climatic shift from historic norms in a given location.For species included in the AquaMaps dataset, we used the thermal preference envelope used to generate species distributions; for species included in IUCN but not AquaMaps, we generated thermal preference envelopes (absolute and preferred minimum and maximum temperatures) in a manner similar to that used to generate envelopes for AquaMaps, using observed mean annual temperature in cells across the species distribution according to IUCN distribution maps: •  $!/ % = (25th percentile -1.5 × interquartile) or absolute minimum mean annual temperature (whichever is lesser) •  $%& % = 75th percentile + 1.5 × interquartile or absolute maximum mean annual temperature (whichever is greater) •  $!/ 0 = 10th percentile of observed variation in mean annual temperature •  $%& 0 = 90th percentile of observed variation in mean annual temperature We modeled physiological thermal stressor intensity  1 for each species based on the local mean annual temperature  ‾ relative to its preferred and absolute thermal range: For each species, mean annual temperature from [75] in each pixel across its distribution was compared to the thermal preferences according to the above formula to generate a species-specific map of thermal stressor intensity.Note that species whose minimum depth preference was deeper than 200 meters (i.e., not epipelagic) were assigned a value of zero for this stressor.

Stressor layers for habitat method
Stressor layers for the Habitat Method of calculating cumulative human impact included all the same uniform exposure stressor layers used for the species and functional entity methods.See S4 Table for details on transformations and rescaling.Species-specific stressor layers (SST rise stressor based on species thermal tolerance, targeted fishing based on species identity, benthic and pelagic bycatch based on position in the water column) could not be included in the same manner as for the species approach, as those calculations depended on species-specific information that cannot be broadly applied at the representative habitat level.However, fisheries pressures (targeted and bycatch) were accounted for by creating layers from the same source, i.e., Watson [33], in the same manner as described in Halpern et al. [7].These new layers were aggregated by method, depth, and scale: commercial pelagic and demersal low bycatch, commercial pelagic high bycatch, commercial demersal destructive, and artisanal/small scale fishing.For all these categories, overall fishing intensity (catch per km 2 ) was normalized by ln(NPP), then rescaled to the 99.9% quantile across all years for that catch category, resulting in stressor scores from 0 to 1. See S4 Table for details on transformations and rescaling.

Habitat maps for habitat method
Habitat maps prepared for Halpern et al. [7], at 934 m resolution in a Mollweide equal area coordinate reference system, were aggregated by a factor of 11 (to identify habitat density at approximately 10 km resolution) then reprojected to match our analysis grid.Corals and seagrass layers were based on habitat maps updated for Berger et al. [79].Kelp and saltmarsh layers were based on maps updated for the Ocean Health Index 2021.

Functional entities
We assigned species to functional entities based on categorical values of four traits (maximum body length, adult mobility, water column position, adult trophic level) that roughly gather species into similar niche space, following Mouillot et al. [26].Out of 512 possible functional entities (8 body

Figure S2 .
Figure S2.Global pattern of species richness based on 21,159 species range maps included in this assessment.

Figure S3 .
Figure S3.Percent of cells falling into each impact quartile category.(A) Species method.(B)Habitat method.Note that a uniform distribution would result in a value of 6.3% in each impact quartile category.

Figure S4 .
Figure S4.Distribution of modeled risk of impact based on vulnerability and exposure of ecosystem-representative habitats to anthropogenic stressors.(A) Mean cumulative impact across all habitats in each cell, summing across all stressors.(B) Mean cumulative impact across all habitats, summing across all climate-related stressors.(C) Mean cumulative impact across all habitats, summing across all non-climate stressors.(D) Bivariate comparison of distributions of climate impacts (orange) vs. non-climate impacts (purple) by quartile within each stressor group.

Figure S5.
Figure S5.Comparison of cumulative, climate, and non-climate stressors by habitat and species methods across 10 km resolution cells within coastal ( !,  ! ) and oceanic ( " ,  " ) portions of 62 representative marine ecological provinces, transformed to percentile ranks relative to global distribution within each impact category.Filled point indicates median value; bars represent interquartile range (IQR, quartile Q1 to Q3); whiskers indicate observations 1.5x IQR below (above) Q1 (Q3) of bar.Outliers omitted from plots for clarity.A) Cumulative impacts by species and habitat cumulative impact methods.B) Climate impacts by species and habitat cumulative impact methods.C) Non-climate impacts by species and habitat cumulative impact methods.

Table S1 :
Species inclusion by phylum and class.N total refers to the number of marine species in each class according to the World Register of Marine Species [28], N included refers to the number of species included in this assessment.

Table S2 :
Vertebrate species inclusion by class and order.N total refers to the number of marine species in each order according to the World Register of Marine Species [28], N included refers to the number of species included in this assessment.

Table S3 :
Traits used to calculate vulnerability to various stressors, according to sensitivity (S), specific adaptive capacity (A), general adaptive capacity (G), and binary exposure modifier (0/1) components.Where marked, different values of that trait change the value of that component.See[19]for full details and methodology.Note that sensitivity to Biomass Removal and SST Rise stressors (noted with *) are set to 1 for all species, but vulnerability to these stressors is moderated by adaptive capacity and exposure modifier.

Table S4 :
Overview of methods and data sources to generate stressor distribution maps.All stressors projected to 10 km Mollweide coordinate reference system.